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Abstract 

We consider evolution of initial disturbances in spatially extended systems 
with autonomous rhythmic activity, such as the heart. We consider the case 
when the activity is stable with respect to very smooth (changing little across 
the medium) disturbances and construct lattice models for description of not- 
so-smooth disturbances, in particular, topological defects; these models are 
modifications of the diffusive XY model. We find that when the activity on 
each lattice site is very rigid in maintaining its form, the topological defects — 
vortices or spirals — nucleate a transition to a disordered, turbulent state. 

PACS numbers: 87.19.Nn, 64.60. Cn 
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I. INTRODUCTION 



Physical mechanisms underlying many cardiac arrhythmias, in particular the transition 
from ventricular tachycardia (VT) to ventricular fibrillation (VF), are not fully understood. 
The ventricular tissue is known, both experimentally and theoretically, to support long- 
living spiral excitations, and it is thought that a breakup of such a spiral could give rise to 
a turbulent, chaotic activity commonly associated with VF. (Spirals are reviewed in books 
IJ.) A considerable effort is now being directed towards understanding of these defect- 
mediated transitions to turbulence within mathematical models of ventricular tissue. The 
currently popular approach (reviewed in Ref. 0) considers a spiral in a patch (or slab) of 
ventricular tissue; the patch is taken in isolation from any pacemaking source. One then 
follows numerically the time evolution of that initial spiral. 

In the real beating heart, however, the ventricles are not isolated from other regions, and 
the heart, viewed as a whole, supports a (more or less) periodic autonomous activity — the 
heartbeat itself. In this case, any defect should be properly viewed as a disturbance of the 
normal heartbeat, rather than a structure in isolated tissue. In this paper we present some 
general results on the evolution of initial disturbances in autonomously active media and 
discuss their possible applications to cardiac arrhythmias. In particular, we identify a simple 
mechanism of defect-induced transition to turbulence in discrete (lattice) systems. We also 
find that the more rigid is the system in maintaining locally the undisturbed form of activity, 
the more easily the transition to turbulence occurs. This observation can potentially identify 
a useful therapeutic target. 

The assumed lattice structure need not (though it may) be related to the mechanical 
structure of the medium. The size of the lattice spacing in our models simply represents 
the smallest spatial scale on which the rhythmic activity can be desynchronized: a region 
smaller than that scale will necessarily fire as one. Discrete models of fibrillation have a 
long history, cf. the 1964 model of Moe et al. ||. (Unlike these authors, though, we do 
not introduce any frozen inhomogeneity in the parameters of the medium, apart from the 
lattice structure itself.) In addition, the importance of a discrete (granular) structure of the 
medium has been emphasized in theoretical studies of defibrillation [Q. 

We introduce an interaction of an excitable region (like the ventricles) with a pace- 
making region using the following simplified (not anatomical) model. We consider a three- 
dimensional (3d) slab of simulated medium whose extent in the z direction is limited by the 
planes z = and z = L z . The properties of the medium change in the z direction: the region 
near z = is spontaneously oscillatory and represents the pacemaking region; the region 
at larger z is merely excitable and represents the ventricular tissue. The z direction will 
be also called longitudinal, and the other two directions, x and y, will be called transverse. 
The medium supports a spontaneous rhythmic activity, in which an infinite train of pulses 
propagates from small to large z. This steady activity is independent of x and y and is 
supposed to model the heart's normal rhythm, in which pulses propagate from the inner 
surface of the ventricles out. 

The goal of our study was to see what happens if at some instant the spontaneous 
rhythmic activity is disturbed in a spatially nonuniform fashion, and then the system is left 
to itself. We approach this question in two steps. First, we consider the case when the initial 
disturbance is very smooth, i.e. almost uniform across the medium; in particular, it captures 
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no topological defects. In this case, we expect that locally the activity rapidly relaxes close 
to its undisturbed form. The state can then be described using a single field r(x, y, z; t), 
which measures the space- and time-dependent delay (or advance) in activity among the 
local regions. This field is a phase variable: it is defined modulo the period T of the steady 
rhythm. For these smooth perturbations, we expect that the dynamics of r at large times 
will be universal: it will be described by an equation whose form (although not the precise 
values of the coefficients) does not depend on the details of electrophysiology or on the 
microstructure of the medium. In particular, this large-time dynamics does not "see" the 
granular structure of the medium. The form of the equation depends on the symmetries of 
the medium at large scales and can be obtained by keeping terms of the lowest order in space 
and time derivatives consistent with the symmetries. For simplicity, we will assume that 
at large scales the properties of the medium are invariant under translations and rotations 
in the x-y plane and that r does not depend on z, i.e. the disturbance is effectively two- 
dimensional (2d). (Recall that z is the direction of propagation of the normal rhythm.) In 
this case, the equation describing the large-time dynamics has the form 

d t 9 = aV 2 2 9 + c(V 2 9) 2 , (1) 

where the phase 9(x,y; t) is related to r via 

6{x,y;t) = 27Tr{x,y 1 t)/T , (2) 

and a and c are coefficients; V2 is the 2d gradient: V2 = (d x , d y ). 
We define a smooth disturbance by the condition 

|V 2 0| < 2tt/L , (3) 

where L = m&x{L x , L y } is the transverse size of the medium. Under this condition, the 
second term in on the right-hand side of ([]]) is much smaller than the first. We keep it 
nonetheless, because it is the leading term that breaks the 9 — > —9 symmetry. As we 
will see, terms breaking this symmetry play an important role in evolution of non-smooth 
disturbances, such as topological defects. So, it is essential to establish that the coefficient 
c is indeed nonzero. For smooth disturbances, though, the second term is unimportant, and 
eq. (JTT) shows that when a > a smooth initial disturbance relaxes back to the uniform 
steady rhythm (9 = const). The relaxation process is ordinary diffusion. 

It is important to provide a derivation of ([I]) from an electrophysiological model. In 
particular, that would supply certain values for the yet unknown coefficients a and c. In 
Sect. 2 we show how 9 (or r) can be defined within such a model. The smaller are gradients 
of 9, the slower it evolves. One might think that, given an electrophysiological model, it 
should be easy to separate away the slow dynamics and obtain, quite generally, a closed 
equation for 9. This task, however, turns out to be far from straightforward, and as of this 
writing we have not been able to obtain a general derivation of ([!]); in Sect. 2 we illustrate 
the nature of the difficulty. 

To establish that the coefficient c is indeed nonzero, we then have resorted to the following 
argument. The simple electrophysiological model that we consider can be driven, by a choice 
of the parameters, to a critical (bifurcation) point, at which the autonomous rhythmic 
activity is extinguished. Near the critical point, the system can be described by a complex 
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Ginzburg-Landau (CGL) model of a complex order parameter whose phase is our time- 
delay field 9. For a smooth, almost uniform, perturbation, the CGL description reduces 
to an equation for 9 alone, and that has the precise form ([I]), with definite values of a 
and c. In particular, we find that a > and c ^ 0. As we move away from the critical 
point and towards the form of activity representative of the normal heartbeat, the CGL 
description ceases to be valid. But as it is difficult to imagine how c would now suddenly 
become identically zero, we assume that the large-time dynamics of 9 is still described by 
(H) with a nonzero c. We also assume that a > 0, so that the uniform state is stable. The 
electrophysiological model that we use is reviewed in Sect. 3, and the CGL description is 
derived in Sect. 4. 

The second step of our program is promoting the above description of smooth pertur- 
bations to a description including not-so-smooth perturbations, in particular, topological 
defects. The latter description will not be universal. The lack of universality means (by 
definition) that the description, and the type of the resulting dynamics, depend on the mi- 
crostructure of the medium. Because no activity can be fine-grained indefinitely, it is natural 
to assume a granular, or lattice, structure. In Sect. 5, we construct lattice models and study 
their dynamics. In Sect. 6 we summarize our results. 

II. DESCRIPTION OF SMOOTH DISTURBANCES 

In this section we want to show how the slow variable 9, or equivalently r, can be defined 
within the context of an electrophysiological model. This variable evolves arbitrarily slow 
in the limit of arbitrarily small gradients; it should not be confused with "slow" recovery 
variables of electrophysiology. Our definition of r works for any medium supporting an 
autonomous periodic activity that is stable with respect to smooth, almost uniform, pertur- 
bations. For definiteness, we consider here an electrophysiological equation of the form 

eg-V 2 g-bV 2 g-F(g,g;z) = 0. (4) 

Overhead dots denote time derivatives, V is the 3d gradient, and e and b are parameters. 
The change in properties of the medium in the z direction is described by the function F, 
which explicitly depends on z. Eq. (|) obtains, for instance, when a medium described by 
the two-variable FitzHugh-Nagumo (FHN) model || is placed in an external static electric 
field (we will show that below). In that case, g is the deviation of the recovery variable of 
the FHN model from the static solution. 

We consider cases when eq. (|j) (or, more precisely, a suitable boundary problem based 
on it) has a periodic in time solution of the form 

g(r,t) = <j>(z,t) . (5) 

For example, this solution may describe a train of pulses propagating in the z direction. 
The periodicity means that <f)(z, t + T) = <ft(z, t) for some period T. Notice that, because of 
the translational invariance of (^) in time, <p(z, t — r) is also a solution of @, for any real r 
(albeit with different initial conditions). We now consider a smooth (in space) perturbation 
of the periodic activity described by (|5|) and assume that a sufficiently smooth perturbation 
relaxes back to the periodic state. After the relaxation has been under way for a while, we 
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expect that deviations of g from are already small — except perhaps in the softest mode, 
associated with the time translation. We thus seek a solution to (§) of the form 

g{r,t) = 4>(z,t-T(r,t))+ X (r,t) , (6) 

where r(r, t) is a slowly changing (on the scale of the period T) function of time: f <C t/T. 
In the limit f — > 0, we should return to the solution merely shifted in time, so in this 
limit x should vanish. Thus, when f is small, \ is also small, although not necessarily slowly 
changing. Because of the periodicity of <fi in time, r(r, t) is a phase variable: at each spatial 
point, it is defined modulo the period T. The condition that the perturbation be smooth 
reduces this ambiguity to a common shift by T in the entire space. 

Note that separation of a perturbation into r and x is n °t completely defined by (^): a 
time-dependent variation in r can be absorbed by a variation in x- This ambiguity can be 
fixed by an additional condition — for instance, by requiring that x is orthogonal to <p with 
respect to a certain inner product. Eq. (§) together with the additional condition will then 
provide a complete definition of the slow variable r. 

Now, let us illustrate the nature of the difficulty that arises when one tries to derive a 
closed equation for r from eq. (f|). We substitute ([]) into (||) and expand the right-hand 
side to the leading order in small quantities — the function x an d the derivatives of r. The 
dependence on x will be contained in an expression of the form M(0)%, where M is a linear 
operator, which acts on x an d depends on <p(z,t — r(r, £)). Because of the translational 
invariance of @ in time, the operator M(<f>) almost annihilates <p(z,t — r(r,t)): 

M(4>)4> « ; (7) 

the approximate equality means an equality up to terms of order of the small quantity d t r. 
If the operator M(<f>) were Hermitean with respect to an inner product of the form 

(Xi,Xs)=f dz I dtw(z, t)xi(z, t)x2(z, t) , (8) 
Jo Jo 

for some fixed weight w(z,t), then taking the inner product of (Q) with <fi would, to the 
leading order, project away x an d produce a closed equation for r. In the case of eq. (f|), 
however, the explicit form of the operator M is 

M(0)X = (td 2 t - V 2 d t - 6V 2 - ^ - ^d t ] X , (9) 

where F is F(<f>, 0; z). This operator is clearly not Hermitean with respect to (|8|) with w = 1, 
and indeed we have not found any weight that would render it Hermitean. Thus, we were 
unable to directly separate the slow dynamics of r from the fast dynamics of x- While it 
seems intuitively clear that the slow dynamics will be described by an equation of the form 
(HI), to establish that the coefficients a and c are indeed both nonzero, we had to resort to 
an indirect method, which we describe below. 
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III. A MODEL OF THE HEARTBEAT 



In this section, we describe in some detail the pacemaking mechanism with which we 
model the heartbeat. This simple model, based on the two- variable FitzHugh-Nagumo 
(FHN) kinetics, will be sufficient for our argument justifying (f|) with nonzero a and c. 

Consider a slab of medium described by a FitzHugh-Nagumo model, 

e^ = V 2 E + f(E)-G, (10) 
r)G 

- = £ - 6G , (ii, 

placed in a static uniform external electric field, such as the field of a parallel capacitor. 
Here E is the transmembrane voltage, G is the recovery variable, e > and b > are 
parameters, and V is the 3d gradient. The direction of the external field is our longitudinal, 
or z, direction, and the slab extends in that direction from z = to z = L z . The boundary 
conditions corresponding to this arrangement are 

dE/dz(0) = dE/dz(L z ) = -T , (12) 

where T is a positive constant — the magnitude of the external field. 

The boundary problem (|10|)-(|T2|) has a static solution, E (z), Gq(z). Deviations from 
the static solution are e(r, t) = E(r,t) — E (z) and g(r,t) = G(r,t) — Go(z). Excluding the 
variable e with the help of (|TT|), we obtain an equation of the form (Q) with 

F{g, g- z) = f(E + bg + g) - f(E ) - g - ebg . (13) 

The explicit dependence of F on z appears through the z dependence of E . 

For a range of T the static solution to (p^0|)-(p^) is unstable, for various choices of f(E), 
with respect to arbitrarily small fluctuations of E and G, and the instability gives rise 
to an unending time-dependent activity |J. This will be our pacemaking mechanism. The 
corresponding linear stability analysis introduces a number of useful definitions, so we briefly 
go over it here. 

Expanding eqs. (|10|)-([l"T|) to the first order in e and g, we obtain 



/ de/dt \_ (\ (VI + ft + f[E (z)) 
{ dg/dt i I 



dz 2 

1 




(14) 



This equation should be supplemented by the boundary conditions 



|(0)-|(I»)-0. (15) 
Consider eigenf unctions ip n (z), n > 0, of the z-dependent operator in (|i~4|), 

/ d 2 



\ dz 2 

with the boundary conditions 



f'[E Q (z)} ijj n {z) = \ n Mz) , (16) 
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(17) 



We assume that the eigenfunctions ip n are real and form a complete orthonormal system on 
L 2 [0,L z ). 

The fields e and g can be expanded in the complete orthonormal system {ip n }'- 

OO 

e(r, t) = J2 M n(r 2 , t)ip n (z) , (18) 

n=0 
oo 

g(r,t) = J2M^,t)Mz); (19) 



n=0 



here r 2 is the two-dimensional coordinate: r 2 = (x, y). Eq. ( |14"D then reduces to the following 
second-order in time linear equation 



A n -V|\. I 



Vn+[b + — l - \vn + - 1 + 6[A„- Vi])v n = 0. (20) 



Eq. fl20|) describes a collection of independent oscillators, one for each value of the integer 
n > and of the 2d wave number k. These oscillators have frequencies squared equal to 
uj\ + bk 2 /e and friction coefficients equal to 7 n + k 2 /e, where 

u 2 n = (l + bX n )/e, (21) 
7n = b + X n /e . (22) 

Assuming that the boundary conditions in the x-y plane allow for the k = mode, we 
conclude that the necessary and sufficient condition for instability is that 

A n < max{-e&, -1/6} (23) 

for at least one of the eigenvalues X n . This condition corresponds to there being a negative 
uj\ or a negative j n , or both. 

The parameter e sets the ratio of time scales characterizing changes in the voltage E 
and in the recovery variable G and is typically small. When e < 1/6 2 , the condition ( |23| ) 
becomes 

\ n < ~eb , (24) 



or equivalently 7„ < 0, where j n is the friction (|22|) . 

The question that we now address is whether the condition (^) is ever satisfied for 
physiologically relevant values of the parameters. We choose e = 0.06, b = 0.7, and f(E) = 
6.75E(E — 0.25)(1 — E), as recommended in Ref. |7j for ventricular tissue with "normal" Na 
and K conductances. The only other parameter (besides J 7 ) that we need to choose is L z , 
the thickness of the slab in the z direction. This represents the thickness of the ventricles 
in our simplified model. We have done numerical simulations with L z = 3.2. For lengths, 
Ref. recommends scaling by a factor of 0.5 cm. A somewhat smaller scaling factor of 
0.2 cm is obtained if we equate the characteristic ("Debye") length £ = 0.57, at which a 
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weak static field gets screened inside the medium, to a realistic value of 1 mm. With either 
scaling, though, L z = 3.2 corresponds to a physical length of order 1 cm. 

To find out if the instability occurs for a given value of J 7 , one can numerically solve the 
boundary problem (|L6| )— (|i~7f) and check the condition (|23|) . Alternatively, one can numeri- 
cally integrate the time-dependent problem fllQ|)-(|l2"D with initial conditions corresponding 
to small fluctuations near the static solution. This second approach also allows one to find 
the form of the time-dependent attractor emerging as the instability is cutoff by nonlinear 
effects, so we have adopted it. For the purposes of this section, it is sufficient to consider 
initial fluctuations that are independent of x and y. Using numerical integrations of ([TDD 



(|T2|) with such initial conditions and with the above values of the parameters, we have found 
that the static solution is stable as long as T < T\ ~ 0.4. The value T\ is the lower critical 
value, at which the static solution first becomes unstable as T is increased. The instability 
persists as long as T\ < T < T 2 but disappears when T reaches the upper critical value 
T 2 « I. 

The form of the time-dependent attractor, which develops from small initial fluctuations 
near the static solution, is qualitatively different for values of T that are close to the upper 
critical field as compared to those elsewhere in the instability window. These two different 
forms correspond to propagating versus nonpropagating activity In the range T\ < 
T < !Fp, where T v is somewhat smaller than JF 2 , the attractor is an unending train of 
pulses propagating in the positive z direction. In our model, this corresponds to the normal 
heartbeat. On the other hand, when T v < T < JF 2 , the development of the instability 
is cut off by nonlinear effects when the deviation from the static solution is too small to 
generate a full-fledged pulse. In this case, the entire attractor lies in the proximity of the 
static solution. As T approaches JF 2 , the activity is extinguished gradually: the closer is T 
to JF 2 , the smaller is the deviation from the static solution. This gradual disappearance of 
activity is reminiscent of a second-order phase transition. 



IV. THE CGL DESCRIPTION 

Near the upper critical field, which from now on we will call the critical point, the fields 
e(r, t) = E(r,t) — E (z) and g(r,t) = G(r,t) — G (z) are small (E and G denote the 
static solution). Expanding the system (|10D (pTT]) in e and g so as to retain the leading 
nonlinearities, we obtain 

c)e 1 1 

<~ = V 2 e + f'(E )e + -f(E o y + -/"'(B„)e 3 - g , (25) 

As it turns out, the effect of the e 2 term is relatively suppressed and is of the same order as 
the effect of the e 3 term. So, we kept both types of terms in eq. (|25|). 
Substituting the expansions (|T8|) ([19D into (|25| ) - (|26|) , we obtain 

^ = u n - bv n ; (28) 
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repeated indices are summed over. Here V2 is the 2d gradient: V2 = (d x ,d y ), X n is the 
eigenvalue of the Schrodinger problem (|16|)- ([T7|) , and a and f3 are defined as 



1 [ L * 

O-nrnm' = -7: / ^/"(£())^fmfc , (29) 



2 Jo 

1 r L * 



Pnmm'm" = ~~ / dzf" {E^ n ^) m ^ m i^) m u . (30) 
O JO 

We stay closely enough to the critical point, so that on that side of it where the static 
solution is unstable there will be only one X n satisfying the instability condition (|23|) . That 
will be Ao- In what follows we only consider cases when e < 1/b 2 . Then, the instability 
condition takes the form 

7o < , (31) 
where 70 = b + Xo/e is the friction coefficient (|2"2"|) for n = 0. The closer the system is to the 



critical point, the smaller is |7o|. We make it small enough, so that the frequency squared 



([2~I1) with n = (and hence with all n > as well) is positive and much larger than 7q: 

^ 2 = l/e-6 2 + 6 7 o»7o • (32) 

The large positive u sets the time scale of rapid oscillations of u n and v n . 

We now want to show that when the system is sufficiently close to the critical point its 
dynamics on time scales of order of and larger than |7o| -1 is described by a 2d complex 
Ginzburg-Landau (CGL) model. The field \I/(r2,t) of this CGL model is defined via the 
expansion 

^o(r 2 , t) = (-r^—e-^ + -^e" 2 ^* + c.c.) + ^^t^ + . . . , (33) 
\ b — iujq b — 2iujq J b 

where the omitted terms are higher harmonics, proportional to the third and higher powers 
of exp(±ia;ot); c.c. means complex conjugate. The coefficients A and C are in principle 
series in but near the critical point \1/ is small, and to the leading order A and Co can 
be regarded as constants, which will be determined later. The definition (|33| ) separates away 
the rapid oscillations with frequency ujq and its multiples and, in this sense, is analogous to 
a transition to the nonrelativistic limit in field theory. 

The CGL description is obtained by substituting (|33| ) into eqs. (|27|) — (|28|) , expanding to 
the third order in and finally retaining only terms that contain exp(±iuot) in powers 
0, 1, and 2. One can verify that terms omitted in (|33|) will not contribute to the resulting 
equation. For instance, terms proportional to exp(±3ia>ot) are of order \1/ 3 ; to convert them 
into terms of lower order in exp(±io;ot) one will need to multiply them by at least one power 
of \I/ or ^t, which will make them of the fourth order in 

The CGL description allows us to consider disturbances of the uniform activity that 
satisfy the conditions 

^ = 0(^ 3 ) , = 0(^ 3 ) . (34) 

These are less restrictive than the smoothness condition (H), which now takes the form 
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|V 2 *|/|*| < 2tt/L . (35) 

In particular, unlike fl34"|), the condition fl35| ) explicitly prohibits topological defects, which 
are centered at zeroes of |^|. Under the more restrictive condition (|35|), the CGL dynamics 
reduces, at sufficiently large times, to dynamics of the phase of \1/ alone. 
To the third order in uq is obtained from (28) and ([33|) as 

Mo (r 2 , t) = Cq&V + \y e - luJot + A V 2 e- 2iuJ0t + - - e- iuJot + c.c. J + . . . , (36) 

where dots again denote higher harmonics. As will be checked a posteriori, v n and u n with 
n > are of order \l/ 2 . 

In this approximation, eqs. (|27|)-(|2~8"D with n = become 

e— = (V 2 , - A )m - f - «ooo«o _ 2a oo^o^ - Axjoo^o > ( 37 ) 
^ = M - bv , (38) 
where v > 0, while for n = > they become 

e-— ^ = -A^ - ^ - c^oo^o > ( 39 ) 

= tij, - bv„ . (40) 

We see that in this approximation the modes with n = v > are damped linear oscillators 
driven by the external force proportional to Uq. For the purpose of calculating u u , it is 
sufficient to take computed to the second order in 

u 2 = 2tf f tf + (Ve -2 ^ * + c.c.) + 0(^ 3 ) . (41) 
Then, the solution for u v at large times is 

u u = A u V 2 e- 2iuJot + Al{^) 2 e 2iUQt + C v ¥$ + 0(^ 3 ) , (42) 



where 



A v = -a v00 (x u - 2ieu + j—^j ) ' ^ 
C v = -2a u00 (X u + 1/by 1 . (44) 



Substituting this expression for u v into eq. Q37D for uq we see that the only effect of the 
modes with n > is a local (in space and time) renormalization of the dynamics of the 
n = mode. 

To complete our derivation of the CGL description, we now turn to eq. (^) and compose 
separate equations for different powers of exp(— iuj^t). The equations for the zeroth and 
second powers give expressions for Cq and A that are of the same form as ([43]) -([44]) but 
with v everywhere replaced by 0. The equation for the first power then gives the CGL 
equation 
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* = DVl^ - ^7o^ - , (45) 



where the complex diffusion coefficient is 

1 

2e V ujq 



D=Ui + *L), (46) 



and the complex coupling constant is 

s = B (- 2 |««»(H + ( A »- 2i -» + ^)") +3 A«») ' < 47 > 

Recall that the condition of instability of the static solution is 7 < 0, and near the critical 
point |7q| is small. 

Spatially uniform activity near the critical point (for 7 < 0) is described by the following 
solution of fl4"5[) : 

= Poexp(-is lP 2 t) , (48) 

where p — ( l^o I / 2s/?) x / 2 ; sr and sj are the real and imaginary parts of s. Of course, this 
solution exists only when sr > 0. For a smooth perturbation of this uniform activity (which, 
in particular, contains no topological defects), we can define the modulus p(r 2 ,t) and the 
phase 9(r 2 , t) via 

#(r 2j t) = p(r 2 , t) exp(-«ipgt + 0(r 2j £)) . (49) 



Substituting this into eq. fl33|) shows that 6 1 measures the phase shifts in periodic activity 
among local regions, so it is precisely the variable that we defined in Sect. 2. As the modulus 
p relaxes close to p ~ p everywhere in the 2d space, eq. fl4"5|) reduces to an equation for the 
phase 9 alone. That equation is of the form ([!]), with a = KeD, and c = —ImD. 



V. CONSTRUCTION OF LATTICE MODELS 

As we move away from the critical point and towards the form of activity that is more 
representative of the normal heartbeat, the CGL description ceases to be valid. Nevertheless, 
we expect that eq. ([!]) will still apply for sufficiently smooth perturbations. That is because 
9 is the only variable that can change arbitrarily slowly (for arbitrarily small gradients) , and 
the two terms on the right-hand side of flU) are the only two terms of the lowest (second) 
order in gradients that are consistent with the symmetries of our model and the assumption 
that 9 does not depend on z. Moreover, we now have a reason to believe that both coefficients 
a and c will be nonzero: we have seen that they were both nonzero near the critical point, 
and it is hard to imagine how either of them would vanish identically when we move away. 
So, we consider eq. (P to be reasonably well justified. 

The next step is to build upon ([!]) to construct models that would apply to not-so-smooth 
perturbations of the normal rhythm, in particular, to those containing topological defects. 
As we consider perturbations of progressively smaller spatial scales, there are two effects 
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that lead to deviations from ([[]). On the one hand, the granular (lattice) structure of the 
medium becomes important; on the other hand, the local form of activity deviates from its 
unperturbed form, so that other variables besides 9 come into play. We have found that 
the resulting dynamics depends crucially on which of these two effects becomes important 
first, i.e. at larger spatial scales. In what follows, we contrast the corresponding two types 
of the dynamics. Finding out which one is realized in a specific medium will require a 
detailed electrophysiological model. The required model will have to include the details of 
the granular structure, so it cannot be a simple continuum model of the type we used to 
justify eq. ([[]). 

First, consider the case when the local activity is very rigid in maintaining its form. That 
means that each grain — or lattice site — still carries on essentially the undisturbed activity, 
so the field 9 remains the only requisite variable. In this case, the dynamics is described 
by a model of classical lattice XY spins. For definiteness, we consider here a model on a 
square lattice, with interactions restricted to the nearest neighbors (NN). (Similar results 
were obtained for a model that includes interactions of next-to-nearest neighbors.) We take 
the model equation in the form 

d t 9i = h~ 2 J2 l a sin ( 6 i ~ °i) + C ( X - cos (% - ■ ( 50 ) 

jgNN(i) 

The index % labels the sites of a 2d square lattice, and h is the lattice spacing. Matching to 
the long-wave limit (P identifies a and c in (ID with those in (0). 

Near the critical point, c/a = —b/ujQ, which is proportional to the small \fe. Away from 
the critical point, however, there is no reason to expect \c/a\ to be small, and we need to 
explore the dynamics of the model for diverse values of this ratio. We assume that a > 
and set a = 1 by a rescaling of time. 

When c = 0, eq. (pOD becomes the usual diffusive XY model. This model has stable 
topological defects — vortices and antivortices. A nonzero c gives these defects a rotation 
(clockwise or counterclockwise, depending on the sign of c), so vortices and antivortices 
become spirals. By numerically integrating (|50D, we have found that for small values of |c| 
these spirals are stable — or at least no instability could be detected during finite times of 
our computer runs. 

As |c| is increased, the spirals become more tightly wound and at a sufficiently large |c| 
they become unstable. Formation of a tightly wound but still stable spiral is illustrated by 
Figs. |], H Fig. m shows an initial state, containing a single vortex, and Fig. ^ shows the 
spiral that develops from that initial state for a = 1 and c = —0.5. The values of 9 at a given 
time are represented as directions of lattice spins, as measured clockwise from 12 noon ||. 
These results were obtained via Euler's explicit time-stepping scheme on a 33 x 33 lattice 
with side length L = 10 and discretized Neumann boundary conditions. For picture clarity, 
only a 22 x 22 square is shown. 

Evolution of an unstable defect is illustrated by Fig. [| This picture was obtained for 
a = 1 and c = —2 on the same lattice and with the same initial condition as Fig. |2|. The 
center of the defect now serves as a nuclei of a new phase, a featureless turbulent state. A 
bubble of the new phase originates at the center of the defect and rapidly grows, eating up 
the "normal" phase, until the new phase occupies the entire volume. As far as we can tell, 
the resulting turbulent state is persistent. Fig. ^| shows the bubble during its growth. This 
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growth is indeed so rapid that the initial vortex does not have time to fully develop into a 
spiral, although some fragments of spiral structure can be seen near the wall of the bubble. 
A patch of the turbulent state is seen inside the bubble, away from the wall. When the 
turbulent state occupies the entire volume, it remains disordered: directions of the spins are 
uncorrelated beyond a few lattice spacings. In addition, spins in the turbulent state rapidly 
change their directions with time. 

Next, we consider a case when the local activity is flexible, i.e. it readily changes its form 
in response to a short-scale perturbation. For instance, we can supply the lattice spins with 
a variable length by making 9 the phase of a complex field $ = |$| exp(i6). This introduces 
an additional degree of freedom associated with |$|. As an illustration, consider $ that 
obeys a complex Ginzburg-Landau (CGL) equation: 

^ = DV 2 $ + r$(l-|$| 2 ) , (51) 

where D — a — ic; for simplicity we take the coupling r to be real: r > 0. We can now 
discretize eq. (ETf) on a 2d square lattice of spacing h and vary the parameter r in relation 
to h~ 2 . At large r, the modulus |$| freezes out at |<&| ~ 1, and we obtain a lattice model of 
9 alone, in the spirit (although not necessarily of the exact form) of eq. flSHp. At small r, the 
natural size of a defect's core will be set by (\D\/r) l l 2 , rather than the lattice spacing, so 
we expect that the discretization will be irrelevant, and the dynamics will approach that of 
the continuum 2d CGL model. This latter model has spiral solutions that are at least core- 
stable in a certain range of its parameters |§. Numerically integrating discretized eq. (|51~1). 
we have found that by varying r, for a fixed c/a, one can interpolate between the unstable 
spirals of a lattice model with fixed-length spins and the stable spirals of the continuum 
CGL model. 



VI. CONCLUSION 

In this paper we tried to implement consistently the idea that a disturbance in the normal 
heartbeat can be viewed as a collection of "clocks", each of which measures the local phase 
of the activity. In conjunction with the view that the heart has a granular (or lattice) 
structure, this idea leads to a description of the heart via lattice models of classical spins. 
Our main results are as follows. 

(i) Assuming that sufficiently smooth (almost uniform across the medium) disturbances 
of the normal rhythm relax back to it, one can write down a universal description of this 
relaxation process. Universality means that the form of the equation is independent of details 
of microscopies. For a simplified model of the heartbeat, and disturbances depending only on 
the transverse (with respect to the direction of pulse propagation) coordinates, the universal 
description is eq. ([[]). Although we have not derived this equation in the general case, we 
have justified it by presenting a derivation near a critical (bifurcation) point. 

(ii) For not-so-smooth disturbances, including topological defects, dynamics begins to 
depend on the assumed lattice structure and the details of electrophysiology. In particular, 
we have found that it depends strongly on how rigid the local activity is in maintaining 
its form. When the activity is very rigid (fixed length spins), the system, for a range of 
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the parameter space, is prone to a defect-induced instability, which leads to a disordered, 
turbulent state. 

We expect that the local rigidity of the medium (in the above sense) will depend on its 
longitudinal size (the thickness of the ventricles) and on the electrophysiological parameters, 
such as Na and K conductances. Since, according to our results, the local rigidity plays 
such an important role in the transition to turbulence (fibrillation), its dependence on the 
parameters may serve to identify useful therapeutic targets. 
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FIG. 1. Field distribution at t = 0. 




FIG. 2. Field distribution at t = 20 in the model (|50T) with a = 1 and c 
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FIG. 3. Field distribution at t = 0.3 in the model (|50|) with a = 1 and c 
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